/*******************************************************************************
Continental_US_Grids.do

This file divides up the contiguous US (i.e. 48 states plus DC, excluding
Alaska and Hawaii) into a grid of M-by-M mile squares, for M=3, 6, 12, 24, 48. 
First, we load in state shapefiles obtained from the Census, drop Alaska and 
Hawaii, and combine the state polygons into one polygon for the contiguous US.
Next, we use a projection to convert latitude-longitude coordinates into miles.
We use the Albers equal-area ellipsoid projection, which is used in Holmes and 
Lee (2010, "Cities as Six-by-Six-Mile Squares") and is the generally recommended
projection for areas that are wide but relatively short, such as the United 
States. After projection, we divide the continential US into M-by-M mile 
squares using the "spgrid" command. Because the resulting shapefile of grid 
cells is also given in projected coordinates, we finally invert the projection
to get back to latitude-longitude coordinates.

Last updated: 3/27/2019
*******************************************************************************/

version 15
cd "C:\Plants_in_Space"
set more off
set type double
set varabbrev off

**************************COMBINING STATE SHAPEFILES****************************
// Load in shapefile for the US, downloaded from 
// https://www.census.gov/cgi-bin/geo/shapefiles/index.php?year=2018&layergroup=States+%28and+equivalent%29. 
// This shapefile is at the state-level.
shp2dta using "Shapefiles\Raw\tl_2018_us_state.shp", ///
	database("Shapefiles\Intermediate\2018_state.dta") ///
	coordinates("Shapefiles\Intermediate\2018_state_coordinates.dta") ///
	replace
// Remove areas from outside the contiguous US: Alaska, Hawaii, American Samoa, 
// Commonwealth of the Northern Mariana Islands, United States Virgin Islands, 
// Guam, and Puerto Rico
use "Shapefiles\Intermediate\2018_state.dta", clear
drop if _ID==41 | _ID==32 | _ID==42 | _ID==37 | _ID==36 | _ID==50 | _ID==35
save "Shapefiles\Intermediate\2018_state.dta", replace
use "Shapefiles\Intermediate\2018_state_coordinates.dta", clear
drop if _ID==41 | _ID==32 | _ID==42 | _ID==37 | _ID==36 | _ID==50 | _ID==35
save "Shapefiles\Intermediate\2018_state_coordinates.dta", replace

// Combine state polygons together into one polygon for the contiguous US
use "Shapefiles\Intermediate\2018_state.dta", clear
mergepoly using "Shapefiles\Intermediate\2018_state_coordinates.dta", ///
	coor("Shapefiles\Intermediate\2018_continental_US_coordinates.dta") replace

***************************PROJECTING COORDINATES*******************************
// Reported distances are in latitude-longitude coordinates. We will convert to
// miles.
use "Shapefiles\Intermediate\2018_continental_US_coordinates.dta", clear
// Use an Albers equal-area ellipsoid projection.	
geo2xy _Y _X, gen(new_Y new_X) projection(albers) // Units here are in meters

// Convert meters to miles--note that Y is latitude, X is longitude
replace new_Y=new_Y*0.000621371
replace new_X=new_X*0.000621371

******
// When we invert the projection to get back to latitude-longitude coordinates 
// later on, we need to keep the parameters calculated from this projection. This
// section of code does exactly that, using the instructions given in Snyder 
// (1987), on pages 292-293: https://pubs.usgs.gov/pp/1395/report.pdf.
preserve
	return list
	gen a=r(a) // semi-major axis of reference ellipsoid
	gen f=r(f) // inverse flattening of reference ellipsoid
	gen lat0=r(lat0) // projection's origin
	gen lat1=r(lat1) // 1st standard parallel
	gen lat2=r(lat2) // 2nd standard parallel
	gen lon0=r(lon0) // Central meridian
	destring a f lat0 lat1 lat2 lon0, replace
	replace a=a*0.000621371 // a is originally given in meters, need to convert to miles
	// See Snyder page (13) for relationship between f and e--we take the inverse
	// of that relationship--also see the documentation for the geo2xy_albers help file
	gen e_inv=sqrt(1-(1-(1/f))^2)
	// From Snyder (101, equation 3-12): calculate q0, q1, and q2
	// NOTE: Stata takes values as being in radians, so need to convert to degrees
	gen q0=(1-e_inv^2)*(sin(lat0*_pi/180)/(1-(e_inv^2)*sin(lat0*_pi/180)^2) - (1/(2*e_inv))*ln((1-e_inv*sin(lat0*_pi/180))/(1+e_inv*sin(lat0*_pi/180))))
	gen q1=(1-e_inv^2)*(sin(lat1*_pi/180)/(1-(e_inv^2)*sin(lat1*_pi/180)^2) - (1/(2*e_inv))*ln((1-e_inv*sin(lat1*_pi/180))/(1+e_inv*sin(lat1*_pi/180))))
	gen q2=(1-e_inv^2)*(sin(lat2*_pi/180)/(1-(e_inv^2)*sin(lat2*_pi/180)^2) - (1/(2*e_inv))*ln((1-e_inv*sin(lat2*_pi/180))/(1+e_inv*sin(lat2*_pi/180))))
	// From Snyder (101, equation 14-15): calculate m1 and m2
	gen m1=(cos(lat1*_pi/180))/(sqrt(1-(e_inv^2)*((sin(lat1*_pi/180))^2)))
	gen m2=(cos(lat2*_pi/180))/(sqrt(1-(e_inv^2)*((sin(lat2*_pi/180))^2)))
	// From Snyder (101, equation 14-14): calculate n
	gen n=(m1^2-m2^2)/(q2-q1)
	// From Snyder (101, equation 14-13): calculate C
	gen C=m1^2+n*q1
	// From Snyder (101, equation 14-12a): calculate rho_0
	gen rho_0=(a*sqrt(C-n*q0))/(n)
	save "Shapefiles\Intermediate\inverting_Albers_projection_parameters.dta", replace
restore
******

drop _Y _X
ren (new_X new_Y) (_X _Y)
keep _ID _X _Y
order _ID _X _Y
save "Shapefiles\Intermediate\2018_albers_coordinates.dta", replace

******************************CREATING GRID*************************************
// We use the spgrid command to divide the contiguous US up into M-by-M mile
// squares.
foreach m of numlist 3 6 12 24 48 {
	spgrid using "Shapefiles\Intermediate\2018_albers_coordinates.dta", resolution(w`m') ///
		shape(square) unit(miles) cells("Shapefiles\Intermediate\US_`m'_mile_square_cells_projected.dta") ///
		points("Shapefiles\Intermediate\US_`m'_mile_square_points_projected.dta") dots replace
} 

*****************INVERTING GRIDS TO LATITUDE-LONGITUDE COORDINATES**************
// Each grid is saved in terms of miles in Cartesian coordinates. We need to 
// convert these back to latitude-longitude coordinates, so that we can save 
// each one. For each value of M, we do this to both files created by spgrid: the
// "cells" file and the "points" file.

foreach m of numlist 3 6 12 24 48 {
	******Cells file
	use "Shapefiles\Intermediate\US_`m'_mile_square_cells_projected.dta", clear
	// Merge to get data on parameters in Albers projection, keep all the constants from the merge
	merge 1:1 _n using "Shapefiles\Intermediate\inverting_Albers_projection_parameters.dta", keepusing(a f lat0 lat1 lat2 lon0 e_inv q0 q1 q2 m1 m2 n C rho_0) 
	drop _merge
	drop if missing(_ID)
	foreach var of varlist a f lat0 lat1 lat2 lon0 e_inv q0 q1 q2 m1 m2 n C rho_0{
		sort `var'
		replace `var'=`var'[_n-1] if missing(`var')
	}
	gen rho=sqrt(_X^2+(rho_0-_Y)^2)
	// From Snyder (102, equation 14-11): calculate theta, need to generate this as a 
	// variable since it's a function of each point's coordinates
	gen theta=(atan(_X/(rho_0-_Y)))*180/(_pi)
	// From Snyder (102, equation 14-21): calculate beta, authalic latitude
	gen q=(C-(rho*n/a)^2)/n
	// From Snyder (101, equation 14-9): calculate lambda, longitude
	gen longitude=lon0+theta/n 
	// From Snyder (101, 3-16): calculate phi, latitude, using iterative process as
	// illustrated in Snyder (294)
	gen trial_phi=(asin(q/2))*180/(_pi)
	gen phi=.
	forv k=1/10{
		replace phi=trial_phi+(((1-(e_inv^2)*((sin(trial_phi*_pi/180)^2)))^2)/(2*cos(trial_phi*_pi/180)))*(q/(1-e_inv^2)-sin(trial_phi*_pi/180)/(1-(e_inv^2)*((sin(trial_phi*_pi/180)^2)))+(1/(2*e_inv))*ln((1-e_inv*sin(trial_phi*_pi/180))/(1+e_inv*sin(trial_phi*_pi/180))))*(180/(_pi)) 
		replace trial_phi=phi
	}	
	ren phi latitude
	keep _ID longitude latitude
	order _ID longitude latitude
	ren longitude _X
	ren latitude _Y
	save "Shapefiles\Final\US_`m'_mile_square_cells_latitude_longitude.dta", replace
	
	******Points file
	use "Shapefiles\Intermediate\US_`m'_mile_square_points_projected.dta", clear
	// Merge to get data on parameters in Albers projection, keep all the constants from the merge
	merge 1:1 _n using "Shapefiles\Intermediate\inverting_Albers_projection_parameters.dta", keepusing(a f lat0 lat1 lat2 lon0 e_inv q0 q1 q2 m1 m2 n C rho_0)
	drop _merge
	drop if missing(spgrid_id)
	foreach var of varlist a f lat0 lat1 lat2 lon0 e_inv q0 q1 q2 m1 m2 n C rho_0{
		sort `var'
		replace `var'=`var'[_n-1] if missing(`var')
	}
	gen rho=sqrt(spgrid_xcoord^2+(rho_0-spgrid_ycoord)^2)
	// From Snyder (102, equation 14-11): calculate theta, need to generate this as a 
	// variable since it's a function of each point's coordinates
	gen theta=(atan(spgrid_xcoord/(rho_0-spgrid_ycoord)))*180/(_pi)
	// From Snyder (102, equation 14-21): calculate beta, authalic latitude
	gen q=(C-(rho*n/a)^2)/n
	// From Snyder (101, equation 14-9): calculate lambda, longitude
	gen longitude=lon0+theta/n 
	// From Snyder (101, 3-16): calculate phi, latitude, using iterative process as
	// illustrated in Snyder (294)
	gen trial_phi=(asin(q/2))*180/(_pi)
	gen phi=.
	forv k=1/10{
		replace phi=trial_phi+(((1-(e_inv^2)*((sin(trial_phi*_pi/180)^2)))^2)/(2*cos(trial_phi*_pi/180)))*(q/(1-e_inv^2)-sin(trial_phi*_pi/180)/(1-(e_inv^2)*((sin(trial_phi*_pi/180)^2)))+(1/(2*e_inv))*ln((1-e_inv*sin(trial_phi*_pi/180))/(1+e_inv*sin(trial_phi*_pi/180))))*(180/(_pi)) 
		replace trial_phi=phi
	}	
	ren phi latitude
	keep spgrid_id spgrid_xdim spgrid_ydim spgrid_status spgrid_mapid longitude latitude
	order spgrid_id spgrid_xdim spgrid_ydim spgrid_status spgrid_mapid longitude latitude
	ren longitude spgrid_xcoord
	ren latitude spgrid_ycoord
	save "Shapefiles\Final\US_`m'_mile_square_points_latitude_longitude.dta", replace
}

